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O ' Abstract. During the exothermic adsorption of molecules at solid surfaces 

, dissipation of the released energy occurs via the excitation of electronic and 

phononic degrees of freedom. For metallic substrates the role of the non- 
adiabatic electronic excitation channel has been controversially discussed, as the 
absence of a band gap could favour an easy coupling to a manifold of electron- 
hole pairs of arbitrarily low energies. We analyse this situation for the highly 
exothermic showcase system of molecular oxygen dissociating at Pd(lOO), using 
, time-dependent perturbation theory applied to first-principles electronic-structure 

G ■ calculations. For a range of different trajectories of impinging O2 molecules 

' we compute largely varying electron-hole pair spectra, which underlines the 

' necessity to consider the high-dimensionality of the surface dynamical process 

(-H ' when assessing the total energy loss into this dissipation channel. Despite the 

high Pd density of states at the Fermi level, the concomitant non-adiabatic energy 
losses nevertheless never exceed about 5 % of the available chemisorption energy. 
While this supports an electronically adiabatic description of the predominant 
heat dissipation into the phononic system, we critically discuss the non-adiabatic 
^■f^ ' excitations in the context of the O2 spin transition during the dissociation process. 
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1. Introduction 



Energy conversion at interfaces is at the centre of the rapidly growing field of basic 
energy science. This concerns desired conversions like solar to chemical energy, but 
also unavoidable by-products like the dissipation of chemical energy into heat. An 
atomistic understanding of the involved elementary processes is in all cases only just 
emerging, but is likely to question established views and macro-scale concepts. With 
respect to the dissipation of heat freed during exothermic surface chemical reactions 
such a prevailing view is that of a rapid equilibration with the local heat bath provided 
by the solid surface. The view is nurtured by the rare event dynamics resulting from 
the typically sizeably activated nature of surface chemical processes: While the actual 
elementary processes themselves take place on a picosecond time scale, times between 
such rare events are orders of magnitude longer. The understanding is then that in 
these long inter-process time spans any released chemical energy is rapidly distributed 
over sufficiently many surface phononic degrees of freedom to warrant a description 
in terms of a mere heat bath with defined local temperature. At the atomic scale 
this equilibration with the surface heat bath leads to an efficient loss of memory 
of the adsorbates about their history on the solid surface between subsequent rare 
events. This motivates the description in terms of a Markovian state dynamics that 
is e.g. underlying all present-day microkinetic formulations in heterogeneous catalysis 
[U [21 [5]. In turn, the freed reaction energy only enters the determination of the 
local heat bath temperature, commonly achieved through a continuum heat balancing 
equation [l|5l|6]. 

In many cases this prevalent frame- 
work seems to allow for an accurate ac- 
count e.g. of catalytic conversions. How- 
ever, particularly for nanostructured sur- 
faces and highly exothermic reactions it 
is presently unclear whether this effective 
description is sufficient. With respect to 
the exothermicity this suspicion comes 
from recalling that freed enthalpies can 
well be of the order of several eVs. This 
is e.g. generally the case for the disso- 
ciative adsorption of molecular oxygen 
at catalytically relevant transition metal 
surfaces, which is the specific surface 
chemical reaction we want to focus on 
here. Several eVs is an enormous amount 
of energy on the scale of phononic de- 
grees of freedom, which calls for effi- 
cient dissipation channels to achieve the 
assumed quasi-instantaneous local equi- 
libration. At metallic substrates elec- 
tronic excitations could hereby poten- 
tially represent an important additional 

energy sink. Unfortunately, no consensus has hitherto been reached concerning the 
role of this additional channel for the gas-surface dynamics and subsequent energy 
dissipation. Among others TuUy and coworkers have frequently stressed that the lack 




Figure 1. Schematic illustra- 
tion of electron-hole pair exci- 
tation during the impingement 
of a gas-phase molecule on a 
metal surface. An electron is 
excited from an occupied state 
below the Fermi level ep to an 
unoccupied state above, result- 
ing in an excited state com- 
monly referred to as electron- 
hole (e-h) pair. 
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of a band gap in metallic systems should in principle allow for electronic excitations, 
namely electron-hole (e-h) pairs, of arbitrarily low energies to easily couple with the 
nuclear motion of a particle impinging from the gas-phase [71 [HI HI [IHl [HI [HI [131 [S] ■ 
This is depicted schematically in figure [H On the other hand, recent comparisons to 
experimental data for hydrogen and nitrogen molecules on metallic surfaces indicate 
that electronically adiabatic descriptions provided within the Born-Oppenheimer ap- 
proximation (BOA) seem to describe the initial interaction dynamics governing the 
dissociation extremely well [HI [HI [JT] . As e.g. recently stressed by Juaristi et al [16] a 
proper high-dimensional description of the nuclear motion appears hereby significantly 
more important than the effect of e-h pair excitations. 

For the particular case of oxygen dissociation the spin-triplet ground-state of 
gas-phase O2 nevertheless brings in another aspect that could require to go beyond 
an adiabatic BOA treatment. With the spin of adsorbed O atoms at metal surfaces 
quenched, a spin transition needs to occur at some stage during adsorption process 
and could then well affect the gas-surface dynamics. The well-known paradigm 
system for this is O2 at Al(lll), where only an explicit account of a hindered 
triplet-singlet transition could reconcile first-principles dynamical simulations with the 
experimentally measured low sticking coefficient for thermal molecules [HI [201 Ell [22] ■ 
This non-adiabatic hindrance was traced back to the inefficiency of both coupling 
mechanisms generally discussed to relax the spin selection rules that suppress reactions 
like oxygen dissociation [191 HQI [211 ISSl [23] : The low mass number of Al leads to a 
small spin-orbit coupling and the low Al density of states (DOS) at the Fermi level 
prevents efficient spin quenching through the preferred tunnelling of minority spin 
electrons between substrate and adsorbate. 

In this context the aim of the present study is to specifically address the electronic 
non-adiabaticity of the O2 dissociation process at Pd(lOO). Apart from the obvious 
catalytic relevance, the specific motivation for the choice of the Pd substrate is 
its extraordinary high DOS at the Fermi level, which is among the highest known 
among the low- index transition metal surfaces [24l [2ll [261 EZ] ■ Much in contrast to 
Al(lll) this and the much higher mass number of Pd should favor an adiabatic spin 
transition, an expectation which we find confirmed by good agreement of a completely 
adiabatically computed initial sticking coefficient with existing experimental data. We 
are going to publish these results in detail elsewhere [28], thereby also addressing the 
most intriguing "different dissociation properties" of cartwheeling and helicoptering 
O2 molecules on Pd(lOO). These have been observed only very recently |29l l30] . 
most strikingly exemplifying the present limits of the aforementioned atomistic 
understanding of an elementary process like oxygen dissociation. Here, however, the 
focus is on the excitation of e-h pairs instead, which according to the general arguments 
of TuUy and coworkers [Zl [9l [HI [Til [El [El E] should be particularly facilitated for 
such a high Fermi-level DOS system. The central objective is thus to assess how much 
of the total amount of ~ 2.6 eV chemisorption energy {vide infra) in this system is 
dissipated into this electronic channel. In addition the computed spin-resolved e-h 
pair spectra will be critically analysed to look for features that relate to the adiabatic 
spin transition. 

2. Theory 

Several approaches of varying accuracy have been suggested to compute electronic 
excitation upon surface adsorption. At the high accuracy end is certainly the " direct" 
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first-principles simulation of e-h pair spectra within time-dependent density-functional 
theory (DFT) and an Ehrenfest dynamics for the nuclei [31] [321 [33] ■ However, 
for the intended semi-quantitative estimate of electronic dissipation, and in view 
of the high computational demands imposed by the O2 at Pd(lOO) system, more 
approximate treatments are appealing. Notwithstanding, the rather rough description 
of the substrate electronic structure in the time-dependent Newns-Anderson model 
developed by Mizielinski et al [34l |35j |36j [37] raises some concerns. The same holds for 
approaches based on electronic friction theory [8l[38l[39l[40], as their results can depend 
sensitively on the way how friction coefficients are calculated [TBI BTl [T51 H31 H5] . 
More severely, already first applications within a forced oscillator model for the 
substrate electrons identify the proper description of spin transitions as an intrinsic 
shortcoming of electronic friction theory [JHlllJ- -f'^'" present purposes we therefore 
opt for an approach introduced recently by Timmer and Kratzer (TK) |48j . which 
relies on perturbation theory applied to a time-dependent DFT framework. As the 
authors have demonstrated convincingly, it does not suffer from any troubles with 
the description of spin, and for the frequently studied reference system of hydrogen 
atoms at Al(lll) very good agreement could be reached with the accurate time- 
dependent DFT/Ehrenfest benchmark data [35]. At the same time the approach is 
computationally very efficient, as it only requires less demanding ground-state DFT 
calculations as input. 

In the following we recapitulate the TK approach. This is done self-contained 
and in considerable detail to highlight several aspects that deserve specific attention 
for the aspired application to the O2 at Pd(lOO) system. Among those is the spin 
transition, which is why also a comparison to electronic friction theory is included. 
Furthermore, several improvements of the original TK approach are described. This 
concerns a numerically more efficient implementation to make the application to 
computationally more demanding systems like the present one tractable. As previous 
work has mainly been focused on purely one-dimensional model trajectories (of mono- 
atomic adsorbates), it also comprises the generalization to curved trajectories (of 
multi-atomic adsorbates). 

2.1. Time- dependent perturbation theory 

The basis for the following derivation is a classical trajectory Q{t) of a particle (a 
single atom or a molecule consisting of N atoms) from the gas phase impinging on the 
metal surface. In general, Q is a 3iV dimensional vector of Cartesian coordinates of 
the atoms, which might be replaced by a suitable reaction coordinate description. It is 
assumed that the particle collides with the surface and is reflected, which is commonly 
referred to as a "full round trip" . Consequently, for very small and very large times t 
the particle is assumed to be inflnitely far away from the surface. This is meant in a 
physical sense such that the particle-surface interaction vanishes. In the philosophy of 
a perturbative approach, the effect of energy loss on the actual trajectory is neglected. 
Accordingly, a rigid surface that does not allow for any phonon excitation is considered, 
and energy dissipation into the electronic degrees of freedom of the surface is assumed 
to be sufficiently small not to significantly change the trajectory itself. 

Thanks to the Runge-Gross theorem [49], the time dependent electronic many- 
body problem can be mapped onto an effective single particle Hamiltonian 

/i^(i) = ic + ^;?D,cffW , (2.1) 
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where to is the kinetic energy contribution and the time-dependent effective potential 
is given by 

vWsit) = vWnit) + «TO,xc(0 + «TD.cxt(i) . (2.2) 

a denotes cohinear spin in two different channels within the framework of spin-DFT for 
collinear spin |50j . v'^^^ describes the electron- ion interaction and explicitly depends 
on time due to the motion of the impinging particle along its trajectory Q{t). The 
dependence of the Hartree potential and exchange-correlation potential on the 
time-dependent density p{r,t) is dropped in order to simplify readability. 

At the "beginning" and the "end" of the full round trip, the (stationary) potential 

^TB.cffW TT; , T:. ^ ^cff:non-int (2-3) 

||Q(t— >±Oo)||2— S-OO 

describes vanishing particle-surface interaction as mentioned before, yielding a (time- 
independent) Hamiltonian 

/i'^W ie+<ff;non-int ^ h^iO) • (2-4) 

So far, no approximations have been made. To obtain a scheme that is more 
computationally tractable than the direct time-dependent DFT approach used e.g. 
by Lindenblatt and coworkers [HIl IMl 133] , the time-dependent effective potential is 
now approximated by its counterparts in a "series of snapshots" of respective separate 
non-time-dependent ground state problems along the considered trajectory Q{t) as 
described by stationary DFT, 

VTB,es(t) « MQ{t)) + <c(Q(i)) + Vo.t{Q{t)) = V'iQit)) + <ff,„o„-i„t , (2-5) 

^ V ' 

=<tf(Q(*)) 

with the interaction part of the approximated effective potential v'^(t) vanishing at 
the boundaries of the considered trajectory 

v^it) '-^^ > 0. (2.6) 

||Q(t^±oo)||2^oo 

Following the arguments of TK [48], the concomitant neglect of deviations from the 
instantaneous ground state charge density can only be justified for systems where the 
densities associated with excited electrons and holes are small compared to the overall 
charge density. From a theoretical point of view, this can, of course, be confirmed a 
posteriori if the calculated electron-hole pair spectra are sufficiently small {vide infra). 
Equations p.4|) and ()2.5|) define an approximation for (|2.ip . 

h''{t)^h^o)+v^{Q{t)) , (2.7) 

which is a quantum mechanical textbook example for a starting point of time- 
dependent perturbation theory: The unperturbed Hamiltonian /i(o) gets acted on 
by a (implicitly) time- (and spin-) dependent perturbation potential v'^{Q{t)). This 
potential might not be small, but should only be non-zero, cf. ()2.6|) . during a short 
time interval, given by the time of "intense interaction" (collision) of the impinging 
particle with the surface. Up to first order, the transition rate pij for an excitation 
within one spin channel a from an occupied state |ef ) into an unoccupied state |ej) 
(i.e. i ^ j), which belong to the part of /i(o) that describes the (clean) surface, is then 
given by 

= ^ (eJb"(QW)kn expQ(4-£ni) • (2.8) 
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States associated with the impinging particle (spanning the other part of the Hilbert 
space which /ij-q) acts on) are not considered here, as the e-h pair excitations to 
be described are supposed to be located within the substrate. For closer contact 
with actual calculations, discrete sets of initial and final substrate states rather 
than continuous band structures are denoted by indices i and j, respectively, in the 
following. In practice, theses indices hence encapsulate fe-point and band indices of 
the Kohn-Sham states of the clean surface. 



2.2. Electron-hole pair spectra 

The excitation spectrum for a complete round trip is obtained from (|2.8|) by integrating 
over time and summing over allowed transitions: 

2 



P,^(fi^) = ^ 



+ 00 



Atpij(t) 



S {hLO - (4 - <)) 



(2.9) 



Here huj are positive and non-zero excitation energies. Furthermore, here and in the 
following, for these kind of sums over transitions appropriate weight factors (depending 
on symmetry) for the fc-point part of the indices i and j are implicitly considered. 
Integrating (|2.9p by parts leads to 

+00 

dtpUt) = 



1 



F ■ — £ 



- exp(-(ej-ef)t 



+ 00 



e — e 





dv'^ 






dt 



exp 



-ie-~sf)t 



(2.10) 



Here, the boundary term vanishes due to the properties of the interaction potential 
v'^ as given by (|2.6I) . Reinserting the result of (|2.10p back into (|2.9p gives 
2 



p,;(fi^) = J2 



£ — £ 



6 {tKo ~ (£j - ef )) , 



(2.11) 



where the "dressed" (transition) matrix elements 



Af,. = / At (ej I VQV^Qit)) I e,r) • Qit) ■ exp ( ^(ej - t 



(2.12) 



are ultimately the key ingredients to be calculated as outlined in section 12.61 below. 

Up to now, only occupations of substrate states corresponding to zero temperature 
have been considered in (|2.8p . (|2.9I) and (|2.11[) . Generalization to finite electronic 
temperatures is straightforward by introducing Fermi occupation factors 

/(e) = . (2.13) 



exp 



£ — £F 



1 
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Not allowing for de-excitations of "thermally smeared" electrons, (|2.1ip then gets 
extended to the following final working expression for the electron-hole pair spectrum: 
2 

5 (ej - el)) [/ {el) - f (ej)] 6 (h^) . (2.14) 



Pc'^i.f^) = E 



e — e 
3 « 



Individual spectra for electrons and holes can be obtained consistently by considering 
respective transitions only relative to the Fermi level, 

2 



nx,ho(M=E 




6 {t^ - {e^^ - s^)) [f{eT)-f{e^,)] 9 (hoj) 



(2.15a) 



S {hLo - (£- - ep)) [/ {ef) - f (ej)] 6 {-hio) , (2.15&) 



where, as before, the excitation energies huj are positive and non-zero. As pointed 
out by TK [48j (|2?T4| . (j2.15ap and (j2.15&P stay mathematically well defined even for 
\ej — Sil — >■ 0. Nevertheless, numerical evaluations have difficulties in these regions, 
which is why they are omitted in plots of calculated spectra below. 



2.3. Comparison to electronic friction theory 

When electronic friction theory [351 IMl US is combined with the forced oscillator 
model |46[ 147] and introducing further approximations to obtain excitation spectra 
[5T1 [52] , the equations directly comparable to (I2.15ap and (|2.156p read as follows when 
adapted to the present notation: 



— / de 



(ep + Sw) — e 



[/(e)-/(e^ + M] ^(M (2.16a) 



-fpr 



FOM; ho 



(huj) 



+ 00 

Trn 



(ep + ^) 



[/(4 + M -/(£)] Oi-huj) , (2.166) 



with 



A^,i,(£) = / dt y/v^iQ{t)) Q{t) exp ( -et 



(2.17) 



based on the electronic friction coefficient 

d<ff'(Q(i)) 



r^%Q{t))^iThY, 



e-{t) 



dQ 



S{eUt) 



mit) 



(2.18) 



For the sake of simplicity and better readability, but without loss of generality, a 
one-dimensional description via a reaction coordinate Q{t) has been employed here. 

Comparing (|2.12p and (|2.17l) the key ingredient to be calculated for both theories 
is very similar: Matrix elements of the type 

d<ff(Q(i)) 



4(0 



dQ 
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with v'^^{Q{t)) as defined in (I2.5p . Obviously, it does not matter for tire derivatives in 
(|2.12p and (|2.18|) that the constant part ^cff- non-inf ^IM^ is not subtracted here. 
In fact, due to the energy differences entering X'J^-^^ in (|2.16qP and (j2.16fep . this can also 
be read as matrix elements for electronic transitions - similar to A^^- . Their dependence 
on the electronic structure is less direct though, i.e. "only" due to the forced oscillator 
model connected to the electronic friction coefficient 77*^. Additionally, it is important 
to note, that instantaneous Eigenstates \e1it)), |ej(i)) are used to obtain the latter. 
Obviously, the present theory can in principle be modified in a straightforward way to 
use instantaneous Eigenstates as well in (|2.12l) . However, this leads to conceptual 
issues with shifts of Eigenvalues in the DFT calculations at different points of a 
trajectory and technical problems with (relative) phase shifts of the instantaneous 
Eigenstates. 

Apart from this, the main difference between the present and electronic friction 
theory comes from the fact that squared moduli of the (respective) matrix elements 
appear before the time integration is carried out in (|2.17p and (|2.18p . As discussed 
extensively by TK [48], this is the reason why ((2l^ . (|2J4)) . (|2.15aP and \2.1hh\ 
remain mathematically well defined even in case of an adiabatic spin transition: The 
spin population assigned to the impinging adsorbate typically (at least for hydrogen 
atoms impinging on different substrates) [311 EH ESI |46l |47j |48] is proportional to 
|Q(t) — QoT^^ around the transition point Qq. This leads to a |Q(t) — QqP^^^ 
singularity of dw^jj ((5(i))/dQ, for which the aforementioned main difference is of crucial 
(mathematical) importance. 

2.4. Dissipated energy 

Finally, the energy dissipated into electron-hole pair excitations, which is the key 
objective here, is obtained as energy-weighted integral over the excitation spectrum. 





Again, a comparable expression is also obtained within electronic friction theory. 



which additionally would allow to "measure" the dissipated energy for only a part 
of the considered trajectory Q{t). However, with friction theory not being (directly) 
applicable here, (|2.19l) should serve well for providing the desired estimate of energy 
dissipated into electron- hole pairs. 

2.5. Realistic first-principles trajectories 

In previous work addressing e-h pair excitation typically strictly one-dimensional 
model trajectories Q{t) along the surface normal have been considered when studying 
adsorption of mono-atomic adsorbates with perpendicular impingement over high 
symmetry sites [3ll[32l|33llMlllIl|48l[5ll[52]. The motion of incoming poly-atomic 
particles (like the diatomic molecule O2) is in contrast less likely to be reducible to 
a single spatial dimension. To keep the input data entirely scalar, an effectively 



(|2T4| : 




(2.19) 




(2.20) 
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one-dimensional description of coordinates and velocities, i.e. a suitable reaction 
coordinate Q, must therefore be employed. A canonical choice is the arc length 
associated with a given trajectory Q(t): 

t 

Q{t) ^ JdT ||AQ(r)||2 . (2.21) 

T 

However, we note that in case of non-negligible changes of relative coordinates within 
the adsorbate, (|2.21|) does not provide immediately obvious physically meaningful 
information about the trajectory of an adsorbate. For a stiff diatomic molecule like 
O2 this might not be much of a concern, unless considering large vibrational excitations 
e.g. during scattering at high kinetic energies. Still, it is interesting to note that (|2.12p 
and all its descendants remain invariant, if trajectories are described by scaled and 
shifted coordinates and consistent velocities 

Q{t) — > Q{t) = Qo + aQ{t) (2.22a) 
Qit) Q{t) = aQ{t) . (2.22&) 

When a = 1 / %/ A^ads for an adsorbate consisting of Aads atoms, the corresponding arc 
length Q{t) then conveniently corresponds to the absolute distance travelled by the 
centre of mass of the adsorbate (or likewise any of its constituents). 

The trajectories employed in the present work have been obtained from molecular 
dynamics simulations on a six-dimensional (Born-Oppenheimer) potential energy 
surface (PES) Veo, describing the molecular degrees of freedom of an isolated 
O2 molecule on a frozen Pd(lOO) surface. This continuous PES representation is 
constructed by a highly accurate interpolation of discrete DFT energies computed 
within the Perdew, Burke and Ernzerhof (PBE) generalized gradient approximation 
to electronic exchange and correlation [53l[54]. Further details and a statistical analysis 
of the resulting dissociation dynamics will be published elsewhere [25 . Briefly, the 
interpolation is based on neural networks, adapted to properly take symmetry into 
account 4000 single-point DFT calculations in supercell geometries, carefully 

sampling molecular positions and orientations over the Pd(lOO) substrate, serve as 
database for the interpolation. The supercell geometry is characterized by (3 x 3) 
supercells, a five layer slab, and a vacuum distance of 15 A. A (4 x 4 x 1) Monkhorst- 
Pack grid is used for fc-point sampling. Electronic states are described by a plane 
wave basis using a cut-off energy of 400 cV together with ultrasoft pseudopotentials 
as implemented in the ICastepI code and supplied within the default ICastepI 
pseudopotential library |56) . 

2.6. Efficient matrix element calculation 

The electronic structure data obtained concomitantly with the single-point energetics 
for the PES then also forms the basis for the evaluation of the matrix elements entering 
(j2.12l) . In several applications of electronic friction theory, the corresponding matrix 
elements given by (j2.17l) are related to the nuclear kinetic energy terms, which are 
neglected within the Born-Oppenheimer approximation. Following ideas from Head- 
Gordon and TuUy [71 18] , the commonly used plane wave implementation by Lorente 
and Persson [S71 [SSI [SI] makes use of first order changes in Kohn-Sham states and 
Eigenvalues to obtain the aforementioned matrix elements. Obviously, this cannot 
work in the present context when dealing with unperturbed states of the substrate. 
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In the original TK work [48^ and in an equivalent strategy previously used within 
electronic friction theory [32] , the matrix elements are calculated directly through the 
derivative of the effective Kohn-Sham potential. Quite in contrast, here we follow 
a different strategy instead, namely to compute M^^ — (ej | v'^^{Q) \ ef) at all points 
Q{tn) for which electronic structure data is available. Afterwards, each of these matrix 
elements is interpolated along Q and the analytical derivative of the interpolation gives 
dM[^{Q)/dQ. Through 



dv'^iQjt)) 
dt 



^ ^M^^iQ)-Qit) , (2.23) 

this then also yields the matrix elements required to determine the e-h spectrum via 
(|2.12p . Following common practise to rely on sufficient decoupling between periodic 
images (SU |32j |33l |48l [57l |58l [5^, only intra-k transitions are considered in this 
procedure, which becomes exact for a periodic overlayer [55]. 

The crucial advantage of this new strategy is that it allows for efficient use 
of existing infrastructure in most plane wave codes, because the derivative of the 
potential does not need to be constructed in a form which can be used to act on 
the states directly. This is particularly useful (if not essential) when dealing with 
non-local and/or ultrasoft pseudopotentials. Notwithstanding, the matrix elements 
calculated for all Q(i„) have to be kept in memory to perform the interpolation and the 
ensuing derivative. This does lead to higher memory demands than for other schemes, 
which by virtue of (memory) parallelisation nevertheless does no longer constitute an 
invincible problem. We have carefully validated and benchmarked our implementation 
strategy by recomputing the H at Al(lll) data reported in the original TK work |48]. 
Quantitative agreement is reached at a minute fraction of the computational cost. 
Similar observations have also been made by Timmer and Kratzer in their more recent 
work [51] after also adapting to our implementation strategy. 

In theory, the loss of orthonormality between substrate states when using ultrasoft 
pseudopotentials [6OI [61] could pose a serious problem for (|2.8I) . hence shaking the 
foundations of the perturbative approach. However, in practice, the aforementioned 
tests for H at Al(lll) have shown that this is not a general concern. Selected 
comparisons of the results for the present system to results obtained with norm- 
conserving pseudopotentials indicate some quantitative dependence of the ultimately 
deduced dissipated energy, yet, without touching on any of the physical conclusions 
discussed below. Much more problematic is the well-known spurious ferromagnetism 
of Pd within semi- local DFT [62j [63] . Test calculations indicate that this leads to a 
dramatic increase of the dissipated energy that can be of the order of several 100 % in 
relative terms. For all results presented below, matrix elements have therefore been 
evaluated using the Kohn-Sham states of a non-spin polarized clean Pd(lOO) surface 
that have been "cloned" into both spin channels of the (spin-polarized) adsorption 
calculations including the O2 molecule. In this respect the fact that only the derivative 
of the potential enters (|2.12l) can be seen as particular (additional) virtue of the 
employed scheme for this specific system: It allows to exploit a cancellation of errors 
in differences (as already in case of the energetics) - quite in contrast to a conventional 
treatment e.g. within time-dependent DFT. 
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Figure 2. Electron- hole pair excitations created by an O2 molecule impinging 
side-on above a hollow site (h-para) as shown in the inset at the bottom, (a) 
PES VgD along the trajectory given by the reaction coordinate Q (neural network 
interpolation = black solid line, DFT input data = black circles), as well as 
projections of the spin density onto the two constituting oxygen atoms (O'^, 
= dotted lines in shades of dark red, sum of O'^ and = light red solid line), 
(b) Evolution of reaction coordinate Q{t) and corresponding velocity Q{t) with 
time t along the trajectory, (c) Separate electron (at positive excitation energies) 
and hole (at negative excitation energies hw) spectra ^'oxol(^) ^^'^ ^Gxho(^) 
according to ||2.156| | and Il2.15al l. (d) Total e-h pair spectrum P^(hu]) given by 
l|2.14|l together with resulting dissipated energies according to l|2.19|l . All spectra 
are for a half round trip with energies hw relative to the Fermi energy. Both 
majority (t, violet) and minority (J,, blue) spin channels are shown. 



3. Results 

3.1. Electron-Hole Pair Spectra and Dissipated Energies 

With the objective to compute e-h pair spectra and to estimate the concomitant 
energy dissipated into electronic excitations we focus our investigation on four distinct 
trajectories, selected to span the space of possible O2 impingement over the Pd(lOO) 
surface. Namely these are two trajectories in which the O2 molecule initially 
approaches with its molecular axis parallel to the surface and oriented along a [100] 
direction (a so-called side-on configuration), and two trajectories in which the O2 
molecule initially approaches with its molecular axis perpendicular to the surface 
(so-called head-on configuration). In order to sample the influence of the lateral 
corrugation of the surface potential one of the two side-on trajectories starts out 
with the O2 center of mass located above a fourfold hollow site, and one trajectory 
with this centre above a twofold bridge site, cf. insets in figures [2] and [H Similarly 
the two head-on trajectories start out above a hollow and above a top site, cf. insets 
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Figure 3. Same as figure [2] but for an O2 molecule impinging side-on above a 
bridge site (b-para) as shown in the inset at the bottom. 




Figure 4. Same as figure [2] but for an O2 molecule impinging head-on above a 
hollow site {h-perp) as shown in the inset at the bottom. 
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Figure 5. Same as figure [2] but for an O2 molecule impinging head-on above a 
top site {t-perp) as shown in the inset at the bottom. 



in figures |3] and [S] Henceforth we will use the intuitive abbreviations h-para, b-para, 
h-perp and t-perp for the four cases, respectively. Each trajectory begins with the O2 
center of mass at a perpendicular distance Zo = 6 A above the surface, where the 
adsorbate-substrate interaction potential, i.e. the PES VgD, has already well decayed 
to zero. As an initial kinetic energy we choose 50 meV to mimic thermal molecules 
in the high-energy Boltzmann tail. In general, non-adiabatic energy losses should 
increase with higher velocities, which is why we expect the deduced dissipated energy 
to represent a useful upper bound for a thermal gas. However, in order to explicitly 
test the dependence on the initial velocity quantitatively we additionally consider a 
h-para trajectory with an increased initial kinetic energy of 400 meV, corresponding 
to an almost tripled initial velocity from Q{0) = 5.5A/fs to Q{0) = 15 A/fs. 

The results for the four trajectories h-para, h-para, h-perp and t-perp with low 
initial kinetic energy are compiled in figures [2] to [5l respectively. In all four cases 
the molecule is reflected, such that a half round trip as assumed in the perturbative 
approach can be consistently defined. With the definition given in (|2.22a|l the reaction 
coordinate Q{t) corresponds for the considered trajectories essentially to the vertical 
distance Z{t) of the O2 center of mass above the surface. As shown in parts (b) 
of figures [2] - [5] the half round trip proceeds in all cases in about 600 fs, in which 
the molecule approaches from its initial height of (5(0) ~ Z{Q) = 6.0 A towards the 
surface. Intuitive in terms of the expected corrugation of the surface potential, the 
point of reflection (defining the end of the considered half round trip) is with Z w 2.6 A 
highest for the t-perp trajectory with the molecule head-on above the top site, and with 
Z « 1.0 A lowest for the h-para trajectory with the molecule side-on above the hollow 
site. In the other two trajectories the reflection occurs at intermediate Z heights. 
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These differences in the point of closest encounter have direct consequences on the 
adiabatic spin transition of the O2 molecule. As illustrated by the projected O2 spin 
density in parts (a) of figures[2]-[S]the molecule approaches the surface only sufhciently 
closely in the h-para trajectory for the spin to be completely quenched at the point of 
reflection. In this case exactly the square-root like decay of the spin density results 
that has been generally observed for hydrogen adsorption before [SU [32l [33l [46l |47l |48] 
and which as described in section [2. 31 would prevent application of electronic friction 
theory. In contrast, in the t-perp trajectory the transition to the singlet state has 
not even really started when the molecule gets reflected. However, the adiabatic spin 
transition of the oxygen molecule depends not only on the penetration depth, but 
also strongly on lateral position and orientation. This is nicely illustrated by the 
trajectories h-para and h-perp in figures [3^ a) andSJa), respectively. In both cases 
the minimum Z-height reached is comparable, but the projected spin density on both 
oxygen atoms remains nearly twice as large for h-para. 

Significant acceleration beyond the initial velocity of Q(0) = 5.5 A/fs only takes 
place in areas of attractive potential at low Z heights just before the molecule hits 
the repulsive Pauli wall. As shown in figure EJb) this is obviously largest for the most 
favourable h-para configuration, which is very close to the main entrance channel for 
dissociation in this system [28]. While the maximum speed reached is in this case 
about five times the initial velocity, essentially no acceleration occurs throughout the 
entire trajectory in the most repulsive t-perp case, cf. figure [5{b). These differences 
have to zeroth order exactly the bearings on the generated e-h pairs and concomitant 
energy dissipation one would expect e.g. from electronic friction theory. As shown 
both separately for electrons and holes in figures [Sfc) -Wis) and combined into one 
e-h pair spectrum in figures (Hd) - [Hd) , the qualitative shape of the spectra is for 
all trajectories about the same. However, the absolute magnitude is much larger for 
the high velocity attractive case h-para, resulting into a total amount of dissipated 
energy that is about an order of magnitude larger than in the other three cases. 
Notwithstanding, the latter three cases illustrate that at closer look also other factors 
like the orientation of the molecule matter: Specifically in the trajectories h-para 
and h-perp in figures [31[a) andlD^a), respectively, the O2 molecule experiences the 
mentioned similarly repulsive interaction potential and reaches comparable maximum 
speed before reflection. Nevertheless, the energy dissipation into e-h pairs is by a 
factor of four less for the latter. This constitutes a nice confirmation of one of the 
key results of the electron friction work of Juaristi et al [16] : The importance of the 
"high" dimensionality of the molecule-substrate-interaction also extends to e-h pair 
excitations, such that a proper assessment of the role of this dissipation channel needs 
necessarily to rely on a representative set of impingement scenarios, not just one model 
trajectory. 

Further insight into the velocity dependence comes from the final h-para 
trajectory, in which the initial kinetic energy has been raised to 400 meV. Despite the 
nearly tripled initial velocity, the corresponding trajectory and spectra (not shown) do 
not change qualitatively with respect to the data shown in figure [2l The total energy 
dissipated into e-h pairs is increased by around a factor of 1.3, resulting in a total 
amount of about 100 meV. Compared to the electronic friction theory results for N2 
on W(llO) [12], the relative increase of the total dissipated energy is thus much less. 
There, for perpendicular incidence, a similar increase of the incidence kinetic energy 
has led on average to an increase in the energy dissipation of more than an order 
of magnitude. At present it is unclear whether these differences are due to the two 
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different systems studied, or whether they indicate potential shortcomings of either 
theory. Direct apphcation of the present theory to the N2 on W(llO) system will 
therefore constitute an interesting next step in our work. 

Among the trajectories considered, the total dissipated energy per half round trip 
is thus largest for exactly the h-para trajectory that is "most interesting" for thermal 
impinging molecules in the sense that it falls very close to the dominant entrance 
channel to dissociation into which almost all molecules get steered [28]. For this h- 
para trajectory the total dissipated energy furthermore seems not very sensitive to 
the specific initial velocity. Considering also the approximate nature of the employed 
perturbative approach we therefore take ^ 100 meV as a reasonable estimate of the 
order of magnitude of the electronic dissipation channel. With the present DFT-PBE 
setup we compute a total adsorption energy released in the adsorption of O2, i.e. in the 
exothermic dissociation into two separate O atoms at the Pd(lOO) surface, of 2.6 eV. 
As such the e-h pair excitation channel does not even take up 5% of this energy, 
which provides an a posteriori justification for the assumption of a weak perturbation 
underlying the present approach. On the one hand, the insignificance of this channel 
for this system is somewhat consistent with the inability to detect chemicurrents 
during the adsorption of O2 on polycrystalline palladium [64]. On the other hand, 
it is remarkably low compared to H impingement over the threefold hollow site of 
Al(lll). The latter leads to a comparable release of chemisorption energy, but the 
computed amount of energy taken up by e-h pairs was one order of magnitude more, i.e. 
1 eV (albeit likely favoured by a penetration of the adsorbate into the first substrate 
layer) [32] . Neither is it thus possible to a priori dismiss electronic excitations in the 
adsorption process at metal surfaces, nor does the difference between the two systems 
follow the picture expected from the much higher Fermi-level DOS in the Pd(lOO) 
system. 

3.2. Spin transition 

At a closer look at the e-h pair spectrum for the h-para trajectory shown in figure |2][d) 
a slight asymmetry between electrons and holes deserves further attention. The fact 
that excitations are stronger in the spin (down) minority channel, and specifically in 
the high-energy wing of the spectrum, is particularly interesting in light of the spin 
transition of the O2 molecule, which, according to figure |2][a), is completed for this 
trajectory before the point of refiection is reached. A similar asymmetry, albeit much 
weaker, also seems to be present for the h-perp trajectory, cf. figure |4l correlating 
with the aforementioned similar degree of spin quenching in this case. Timmer and 
Kratzer have obtained nothing comparable in the e-h pair spectra of H on Al(lll) 
[48] . In contrast, for the same system Lindenblatt and Pehlke [Slj did determine a 
similar abundance of excited electrons and holes in the spin majority and minority 
channel, respectively. However, since the time-dependent DFT based dynamics is 
intrinsically non-adiabatic and the spectra are obtained a posteriori by relaxation 
back to the electronic ground state, this is not surprising. It can be easily understood 
by propagation on excited PESs during the non-instantaneous spin transition. In 
contrast, the present description is adiabatic with respect to the spin transition, and 
e-h pair excitations are evaluated based on unperturbed substrate states. 

In order to put the observed spin asymmetry into perspective with the spin 
transition, a detailed analysis of the O2 molecular states interacting with the Pd(lOO) 
surface along the trajectory Q{t) can be illustrative. For this we evaluate the density 
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Figure 6. DOS including projections onto high-lying O2 molecular orbitals at the 
beginning {t = 0) of the h-para trajectory, when there is still negligible interaction 
with the Pd(lOO) substrate. The upper panel shows the majority spin channel 
(spin t), the lower panel shows the minority spin channel (spin J,). Note that 
at this large molecule-surface distance the O2 molecular states 2a, 2it and 2tt* 
are still rather sharp and are here only broadened by a Gaussian with a width of 
0.3 eV for better visualization. The DOS of the clean Pd(lOO) surface (downscaled 
by a factor of 50) is shown in black, with the occupation of states below the Fermi 
level ep represented through the grey filling. 



of states projected onto a Kohn-Sham state \(pQ (Q)) of spin a of the free O2 molecule 
(molecular PDOS) ^65. 

P^^oA^-^Q)-T.\(^oMKiQ))fHe-en) ■ (3.1) 

n 

It is important to emphasize that here the Kohn-Sham states of the interacting system 
of substrate and molecule |e^((5)) are used for this analysis. The most relevant state 
with respect to overlap with the Pd valence d-band is specifically the degenerate 2tt* 
state [TO] , which in the free triplet O2 molecule is occupied by two electrons in the spin 
majority channel and unoccupied in the spin minority channel. Upon approaching the 
surface, the interaction with the substrate potential starts to broaden the discrete 
molecular levels. Their centre of energy along the trajectory Q{t) is each given by 

/oo 
d££p^,.(£;Q) . (3.2) 
-00 

As shown in figure|6l at the beginning of the calculated trajectory this centre is for the 
majority 2^*"^ level located about 2 eV below the Fermi-level, while it is at about 0.5 eV 
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Figure 7. (a) Evolution of the center of the molecular orbital (e2-^,(Q)) 

along the h-para trajectory, (b) Corresponding number of states projected onto 
this orbital in an energy interval 0.4 eV around the Fermi-level, i.e. around the 
dominant energy range of e-h pair excitations highlighted by the green area in 
panel (a). 

above the Fermi-level for the minority 2tt*^ level. If from inspection of figure El^d) we 
roughly take the dominant fraction of e-h pair excitations to spread over a range of 
0.4 eV around the Fermi- level, we see that it is primarily the spin minority 2tt*^ level 
that falls close to this range. To quantify this as the width of the 27r* level increases 
along the trajectory, we evaluate the number of projected states that lie within this 
relevant energy range as 

/•eF+0.4eV 

N^.,{Q)= / de p^,^,{e;Q) . (3.3) 

JeF-0.4cV 

Figure [7] displays this quantity together with the evolution of the centre of the 2tt* 
level in both spin channels along the h-para trajectory. A much larger number 
of corresponding projected states is clearly visible for the spin minority channel 
specifically in the range of vertical heights ~ 2 — 5A that is most relevant for the 
spin transition. Qualitatively the same result is also obtained when reducing the 
considered energy range around the Fermi-level from 0.4 eV to lower values and thereby 
concentrating on the more and more dominant part of the obtained e-h pair spectrum, 
cf. figure [21 

Altogether we therefore conclude that the observed stronger e-h pair excitations 
in the minority spin channel coincide with a larger amount of substrate states in the 
corresponding energy range, which overlap with the 2Tr*^ orbital of the O2 molecule. 
A tempting interpretation is therefore to see the excess excitations as a reflection of 
the preferred tunneling of excited minority spin electrons that leads to the quenching 
of the oxygen spin-density in the here assumed adiabatic spin transition picture. If the 
electronic structure during the "real" spin transition does not behave in a radically 
different fashion, this connection between the different energetic position of majority 
and minority O2 27r* level and corresponding different overlap with the energy range 
spanned by e-h pair excitations would prevail. In that case, this then provides an 
intuitive mechanism how electron tunneling from the (here forcibly) non-spin polarized 
substrate electronic structure compensates the adsorbate spin. 
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4. Conclusions 

In conclusion we have advanced and employed a first-principles perturbative approach 
to study the energy dissipation into e-h pair excitations during the adsorption of 
O2 at Pd(lOO). Despite the unusually large Pd density of states at the Fermi- 
level, which in principle could render this system a prime candidate for electronic 
non-adiabaticity, we obtain energy losses into this electronic channel that are at 
maximum of the order of 5 % of the total released adsorption energy. In detail, 
this loss and the underlying electron-hole pair spectra differ considerably for different 
O2 impingement configurations, i.e. they depend on the molecular orientation with 
respect to the surface and on the lateral position of impact. By far the largest 
amount of excitations is obtained for a side-on approach above the fourfold hollow 
site. This is also the statistically most relevant configuration as it corresponds closely 
to the entrance channel for dissociation into which the dominant fraction of impinging 
molecules gets steered. While these variations confirm the necessity to account for the 
"high" dimensionality of the molecule-substrate- interaction by considering appropriate 
statistical averages over ensembles of trajectories, we find only a weak dependence of 
the total dissipated energy into e-h pairs on the initial O2 velocity. 

All in all we therefore derive from our calculations a value of ^ 100 meV as a 
reasonable estimate of the order of magnitude for the electronic dissipation channel in 
this system. Even when cautiously considering that the approximations underlying the 
employed perturbative approach lead to a similar underestimation of this quantity as 
was identified for H at Al(lll) (4^, e-h pair excitation thus still plays an insignificant 
role for the dissipation of the large amount of 2.6 eV adsorption energy released during 
dissociative adsorption. This assessment is very much in line with the trends observed 
in other studies beyond single atoms at metal surfaces, no matter if the impinging 
diatomic molecule carries a permanent dipole moment (HCl on Al(lll) [66, ) or not 
(H2 on Cu(llO), N2 on W(llO) 16, ), as well as with the hitherto unsuccessful attempts 
to detect chemicurrcnts in experiments over polycrystalline palladium I64j . 




Notwithstanding, two aspects underlying the present state-of-the-art treatment 
need to be emphasized that could in principle radically question our conclusions on 
the role of e-h pair excitations. First, only non-reactive trajectories at a frozen surface 
are considered. While one trajectory is very close to the actual entrance channel for 
dissociation, this still can not capture possibly important non-adiabatic effects during 
the actual dissociation at a mobile substrate. Second, an adiabatic spin-transition 
from the triplet ground-state of the free O2 molecule to the adsorbed singlet state 
is assumed. While a definite answer to the first point has yet to be found, there 
are good arguments that support the latter assumption. The high Pd density of 
states at the Fermi-level as well as the high Pd mass number favour both coupling 
mechanisms generally discussed to relax the spin selection rules that suppress reactions 
like oxygen dissociation, namely spin-orbit coupling and electron tunnelling between 
substrate and adsorbate. With respect to the latter it is noteworthy that an intriguing 
asymmetry of electron-hole pair excitations in the majority and minority spin channel 
is obtained in the present framework. The excess minority spin excitations span 
thereby an energy range around the Fermi-level in which there is a larger overlap of Pd 
substrate states with the minority 2ti*^ orbital of the O2 molecule. In the here implied 
adiabatic spin transition through efficient electron tunneling, this nicely illustrates the 
mechanism with which the non-spin polarized substrate electronic structure quenches 
the adsorbate spin-density. 
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